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Abstract. We develop a semi-analytical method to reconstruct the lensing potential in quadruply imaged gravita- 
tional lens systems. Assuming that the potential belongs to a broad class of boxy non - elliptical models, we show 
how it is possible to write down a system of equations which can be numerically solved to recover the potential 
parameters directly from image positions and using physical constraints. We also describe a code developed to 
search for solutions of the system previously found and test it on simulated cases. Finally, we apply the method 
to the quadruple lens PG1115+080 which allows us to get also an estimate of the Hubble constant Hq from the 
measured time delay as Hq = 56 jlii km s _1 Mpc" 1 . 
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1. Introduction 

Since the very beginning with the discov- 
ery of the first double lens QSO0957 +561 
(Walsh, Carswell & Weynmann, 1979), gravitational 
lensing has found many cosmological applications (see, 
e.g, Narayan & Bartelmann, 199S and references therein) . 
Measurements of time delays in gravitational lenses can 
be used to give constraints to the value of the Hubble 
constant Hq (Refsdal, 1964), avoiding the usual calibra- 
tion problems that plague local distance estimators (see, 
e.g., Frccdman, 2000). In a gravitational lensing system, 
the ray trajectories of the multiple images have different 
paths and pass through different parts of the gravitational 
potential, so that the light travel time is different for each 
image. In cosmology, it is easy to show that the light 
travel time is inversely proportional to Hq, so that by 
the observed time delays between images it is possible to 
deduce the value of the Hubble constant. However, being 
the distance D = cz/Hq (for z < 1), the gravitational 
lensing potential has a key role in the determination of 
D. Unfortunately, application of this technique has been 
slowed by the difficulties of determining time delays (see, 
e.g., the remarks in Schcchtcr, 2000) and of finding a 
unique model of the gravitational lensing potential since 



there are often sequences of lens models which fit the 
same image positions. 

Up to now, about fifty multiply imaged quasars have 



been discovered (see Kochanek et al., 2000 for a detailed 
list) with fifteen quadruple systems. In all cases, the lens- 
ing system consists of a quasar acting as source and a 
galaxy (or a group of galaxies) acting as lens. These sys- 
tems furnish a good tool to probe the potential of the 
lens galaxies. The advantage with respect to the tradi- 
tional stellar dynamics method is that, by the image po- 
sitions, we can measure the shape and the mass of the 
dark halo well beyond the visible edge of a faraway lens 
galaxy. In this sense, quadruply imaged sytems are bet- 
ter constrained than double quasars since the lens models 
have to fit more image positions and the time delays be- 
tween any pair of images. On the other hand, the possible 
constraints from flux ratios have to be used with more 
caution due to the possible systematic influence of mi- 
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crolensing effects (Chang & Refsdal, 1979). In this paper 
we concentrate only on quadruply imaged systems. 

Starting from the image positions, a significant 
amount of numerical computation is usually required 
to derive the parame ters of lensing potential (see, e.g., 
Schneider et al., 1992; ). The de generacy of the results 
is often not fully explored since a large parameter 
space should be analyzed using a huge amount of com- 
puter time. Previously, several authors have often re- 
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strictcd their studies to isothermal or power - law spher- 
ical models (Evans & Wilkinson, 1998), elliptical models 
(Witt & Mao, 1997) and other simple models (Kassiola & 
Kovner, 1993, 1995) with or without external shear. The 
fully general non-parametric method, e.g. the pixellated 
lens method (Williams & Saha, 2000), is very powerful to 



derive the range of degeneracy in the lens models, but it 
involves significant amount of numerical computation and 
does not provide a clear insight to the relations between 
the characteristic parameters of the lens. For these rea- 
sons, it is still desiderable to find analytical, yet general 
potentials, which allow a quick exploration of the param- 
eter space. First steps in this direction are the analytical 
studies of elliptical potentials, with and without external 
shears due to Witt (1996) and Witt & Mao (1998, 2000). 
Witt, Mao & Keeton (2000) have shown that the esti- 
mate of the Hubble constant from the time delay does not 
depend on the angular shape of the potential provided 
that the model is isothermal and with no external shear. 
Beside, they have also shown how it is possible to ana- 
litically estimate the error in the estimate of Hq due to 
assuming incorrectly an isothermal model with no exter- 
nal shear. A further step has been done by Zhao & Pronk 
(2001) who studied a broad class of analytical models with 
non - elliptical and semi - power - law radial profiles. These 
authors have shown how to reconstruct the lens shape and 
the radial profile parameters directly from image posi- 
tions with a semi-analytical method, but they still have 
to assign a priori the "boxiness" parameter so that their 
method still needs some improvements. The aim of our 
work is to develop a semi-analytical method to recover the 
lensing potential parameters and also the source positions 
using only image positions and some constraints coming 
from physical considerations. As a first step, we consider a 
wide class of lensing potentials similar to the one studied 
by Zhao & Pronk (2000) and neglect the external shear, 
which will be considered elsewhere. 

The outline of this paper is as follows. In Sect. 2 we 
show how it is possible to write down a set of equations 
to recover the source coordinates and all the lensing po- 
tential parameters for a broad class of boxy non elliptical 
potentials using as input only the image postions. Solving 
this system requires numerical methods and that is why 
we have developped a code which we describe in Sect. 3 
where we also test the method itself on simulated cases. 
Sect. 4 is devoted to the application to a real quadruple 
system, PG1115+080 : we reconstruct the lensing poten- 
tial and find the source coordinates and then use these re- 
sults to estimate the Hubble constant from the measured 
time delay between images B and C . Finally, discussion 
of the results and future improvements and applications 
of the method are presented in Sect. 5. 

2. Lensing potential parameters from image 
positions 

Let us choose a rectangular coordinate system (x, y) with 
origin on the lens galaxy centre and axes pointing towards 



West and North respectively. Let (r, 9) be the polar coor- 
dinates being 9 the position angle measured counterclock- 
wise from North; then we have the following coordinates 
transformation : 



x = r sin ( 



y = r cos ( 



The time delay of a light ray deflected by the galaxy 
lensing effect is given by : 



At(r,0) = /rVnoo x 

1 o 



rr, cos ( 



(1) 



being (r, 9) the image position, (r s ,9 s ) the unknown 
source position and ip(r, 9) the lensing potential. In Eq.(^), 
h is the Hubble constant Hg in units of 100 km/s/Mpc 
(i.e. Ho = 100/ikm/s/Mpc), while tioo is a typical time 
delay estimated for a given set of cosmological parameters 
(fi m , fU, ^fc) assuming H n = lOOkm/s/Mpc. We have: 



8irGp cr 
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3#o 2 
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which are the density parameters relative to matter, 
cosmological constant and spatial curvature respectively. 
The typical time delay is defined as : 



(D lDos\ (1 + z l ) 
Tl0 ° = l~^J 



(2) 



with the usual meaning for the angular diameter dis- 
tances Dol,Dos, Dls', z l is the redshift of the lens. 

According to the Fermat principle, the images lie at the 
minima of At, so that the lens equations may be simply 
obtained minimizing At. We get : 
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At = 



At = 



r — r, cos (o — i 



r.sm w - u s = 



dr 

1 dtp 
~r~d9 ' 



(3) 
(4) 



The lensing potential ip{ r , &) is usually made out by 
the sum of two terms; the first one is connected to the lens 
galaxy surface mass density E(r, 9) through the relation : 



V 2 ^ = 2^ 



= 2k 



(5) 



being £ crti = (c 2 /4irG)(D s I D lD ls ) and k the di- 
mensionless surface mass density. 

The second contribution comes from the so - called ex- 
ternal shear which we neglect in this analysis^. In studying 
quadruple systems, one may follow two different, but simi- 
lar, ways. On the one hand, one may give E(r, 9) and then 
solve Eq.(||) for ip(r, 9). This ensures that the lens galaxy 
model has a physical meaning, but it is not always possible 

We will return on this topic in the conclusions and further 
again in a following paper. 
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to find analytical solutions of Eq. (g) . This is a severe prob- 
lem since one does not know a priori the exact values of 
the model parameters and so cannot integrate numerically 
Eq.(|J). On the other hand, one may directly give ip(r, 0), 
determine the lensing potential parameters and then solve 
Eq.(||) to get E(r, 0) (assuming that one also knows the 
redshift of lens and images to evaluate S cr jt). At the end, 
one has to check if the derived model is a physical one. We 
follow this second approach since we are mainly interested 
to reconstruct the lensing potential directly from observed 
image positions using a semi-analytical approach. To this 
end we have first to choose an expression (as general as 
possible) for i/j(r, 0). We limit our attention to the class of 
potentials 



where we have defined 



*P(r,9) 



[1 - 8 cos (20 - 20 p )f = r a T(9) 



(6) 



with the constraint 8 < 1. It is evident the factorization 
of the models which are physically well motivated for the 
following reasons. In Eq.(^|), 9 P is the position angle of the 
main axis of the lensing potential which may be also mis- 
aligned with respect to the main axis of the lens model due 
to the possible effect of the external shear or to the influ- 
ence of other galaxies of the same group. For these reasons, 
it is better to not fix 9 p equal to the value observed for the 
lens galaxy so that it should be treated as an unknown pa- 
rameter of the lensing potential. In Eq. (^j) , (3 is a boxiness 
parameter and it is defined such that T(0) reduces to th e 
usual elliptical form when 0=1/2 ( |Witt fc Mao, 1997|) , 



whilst it gives the simple models of Kassiola & Kovner 
(1995) for a — (3 — 1. Redefining in the correct way the 
parameters, Eq.(||) may describe also the class of models 
studied by Zhao & Pronk (2001) to model the quadru- 
ple system PG1115 + 080. The same authors suggest that 
|/3 1 < 1/2 in order to have physical mass -radius relation 
and show that the flattening of the surface mass density 
of the model obtained solving Eq. (||) for the potential @ 
is given by : 



q K = 



1-8 



1-/3 



4/3 



1)8 



1 



4/3 



1)8 



(7) 



We have checked that for many values of a and (3, in 
the physically interesting range, q K < 1 (i.e. the galaxy 
is oblate or spherical) only if 8 < so that our choice to 
consider only models with 8 < 1 does not exclude physi- 
cally interesting cases. It is also easy to see that the surface 
mass density £(r, 9) scales with the radius as r a ~ 2 , so that 
we may impose the constraint a < 2 in order that E(r, 0) 
is a monotonically decreasing function of the radius. 

Let us now intsert Eq.(^) into Eqs.(||) and (||); the lens 
equations are thus : 



r — r, cos ( 



r s sin (U — i 



,) =ar a - 1 T(0) 



(10 



= r a - 1 h(6)H0) , 



(8) 
(9) 



h(0) = 



28(3 sin (20 - 20 p ) 
l-Scos(20 - 20 p ) 



(10) 



Having four images^], in principle we may try to 
solve the eight equations obtained inserting image posi- 
tions in Eqs.(||) and (g) to determine the potential pa- 
rameters (a, [3, S,0 P ) and the source coordinates (r s ,0 s ). 
Unfortunately this is not possible neither analytically nor 
numerically since, in these equations, a and (3 enter as 
exponent and this causes the system to be nonlinear and 
trascendent. However we will show, in the following, how 
it is possible to algebrically manipulate these equations 
in such a way to get a more handable and numerically 
solvable system. 

As a first step, let us insert the coordinates (ri,0i) of 
image i into Eqs.(||) and (^) and divide side by side; one 
gets 



r s cos ( 







*[l-5cos (2di-26 p ) 
28(3 sin (20, - 20 p ) 



(11) 



This relation does not hold if sin (0, — S ) = and/or 
sin (26 i — 26 p ) — 0. Exploiting the lens equations (||) and 
(^|), it is easy to show that these divergences are very un- 
likely. Suppose that sin (26i — 20 p ) — 0; in this case Eq. (j9|) 
reduces to 

r s sin(0 4 -9 S ) =0 

which is satisfied only if r s = or if sin (6^ — 9 S ) = 0. 
The first solution (r s — 0) means that the source and 
the galaxy are perfectly aligned which is very unlikely 
so that we may discard this case. The other solution 
(sin (0i — S ) = 0) is also very unlikely since it means that 
the source is aligned with one of the images. For these rea- 
sons we may assume that Eq.(ll) holds for all images and 
never diverges. The same equation may be also written for 
images j, k, I; dividing side by side the one for image i by 
the ones for the other images, one gets the following three 
equations : 



r, cos ( 



sin ( 



.) 



1 -<5cos(20i 



-6.) sm(0 z -0 s ) 
20 p ) sin (26 j - 20 p ) 



1 - 8 cos (26 3 - 26 p ) sin (20 l - 20 p ) 

Tj -r s cos (6j -6 S ) sin (6 k - S ) _ 
r k - r s cos (6k - 6 8 ) sin(0i-0 s ) 

1 - £cos (26j - 20p) sin (20 k - 20 p ) 
1 - <5 cos (26 k - 20 p ) sin (20* - 26 p ) 

Ti - r s cos (6i - S ) sin (0/ - S ) 



(12) 



(13) 



r x - r s cos (0/ - S ) sin (0 t - S ) 



2 In the following we add to image coordinates a lower index 
running on i, j, k, I to distinguish among them. 
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1 - 5 cos (20 { - 29 p ) sin (20/ - 26 p ) 
1 - 5 cos (26i - 20 p ) sin (29, - 26 p ) 



(14) 



A,/<5 + fig 
PiS +po 



(19) 
(20) 



Note that a and (3 never appear in Eqs.fll2j), (|13|), (|14J) 
which are thus three independent equations for the four 
unknowns: (r S7 6 s ,S 7 6 p ). The system is not closed (i.e., 

the number of equations is not equal to the number of where the coefficients depend only on the two angles 
variables) and so cannot be solved unless we add a further (*« ' 6 p) and are defined as ^Wov,^ : 
independent relation not involving a and (3. To this aim, 
let us consider again the lens equations for image i; solving 
the second one with respect to r^ x T (6{) and substituting 
into the first one, it turns out 



A 



. , - Tj sin (6i - 6s) sin (26 3 - 29 p ) cos (20; - 29 p ) 



n - r s cos (Oi - u s ) = 

a[l-5cos (29, - 29 p )\ 
26/3sm(28i-28 p ) 



>) 



Writing down the same relation for the image j and 
subtracting from the one for image i, one gets : 

(n - rj) - r s [cos (6i -6s) - cos (9j - 6 S )] = 

-9 S )[1 - <5cos(20, - 20 p )] 



r s a | sm 

26(3 



sin (20; - 29 p ) 

sin (6j - 9 S )[1 - 5 cos (29 j - 20 p )] 
sin (29 j - 29 p ) 



(15) 



-n sin (6j - 6 S ) sin (29, - 29 p ) cos (20., - 20 p ) , 
Hij = r.i sin (6j - 6 S ) sin (20, - 20 p ) 

-Tj sin (Oi ~ 6 S ) sin (20^ - 20 p ) , 
Vij = cos (6j — 6 S ) sin (0j — 6 S ) sin (20j 
- cos (6i - 6 S ) sin (6j - 6 S ) sin (20, 
cos (0, - S ) sin (0j - 6 S ) sin (20 
cos (9j — 9 S ) sin (0; — 9 S ) sin (20j 
(r, - r 3 )a k i sin (20, - 20 p ) sin20j 
(rk - n)(Tij sin (20 fc - 20 p ) sin 20; 



20 p ) cos (20, - 20 p ) 
20 p ) cos (26 j - 20 p ) , 



pi 



- 20 



v) 



Po = (n - r 3 )T kl sin (20, - 20 p ) sin 20, 



The same expression may be written also using images 
k and I instead of i and j; dividing side by side these 
relations, one finally finds out the following equation: 



20 P ) 

- 20 p 

- 20 p 

- 20 p 
20„ 



(n - rj) - r s [cos (9i - 6 S ) - 



COS t/j — i 



(rk - n) - r s [cos (9 k - 6 S ) - 
sin(0;-0 s )[l-<5cos (20, 



cos(0, -0 S )] 
- 20J1 





sin (20* - 20p) 




sin (6j 


-0 s )[l-£cos (29 j 


- 20p)] 




sin (29 j - 20p) 




sin (9 k - 


-0 s )[l-£cos(20 fc 


- 20 P )] 




sin (20 fc - 20p) 




sin (0; 


- S )[1 - 8 COS (20; 


-20p)] 


sin (20/ - 20 p ) 



-(rk - ri)nj sin (20 fc - 20 p ) sin 20; 
having posed for sake of shortness : 
ay = sin (6j - 6 S ) sin (20, - 20 p ) cos (29 j - 20 p ) 

- sin (0, - 6s) sin (20^ - 20 p ) cos (20, - 20 p ) , 

Tij — sin 6i — 6 S sin (20j — 20 p ) — sin 6j — 6 S sin (20; — 20 p ) 

It is now very simple to reduce the system of Eqs.(|l7|), 
(|l8|), (19) and ( p0| ) to a system of three equations in only 
three unkwons; equating side by side those equations one 
finally gets the following closed systems in (0 S , 5, 6 P ) : 



aiS 2 + bid + ci 
( 16 ) a 2 S 2 + b 2 S + c 2 



(21) 
(22) 
(23) 



This equation does not contain neither a nor (3 and a 3 5 + b 3 5 + c 3 =0 
thus can be added to Eqs.(pl), (O) and dTi) to get a 

closed system of four independent equations in the four with the coefficients depending only on (0 S , 6 p ) and defined 



variables (r s ,9 S ,8,9 P ). 

This system may be rewritten in a more convenient 
way. Actually, solving each of the Eqs.(|l2|), (|l3|), ( pT[ ) and 
( p^| ) with respect to r s we simply get : 

Ay 6 + fiij 



v l j8 + rjij 

A;fc<5 + [Hk 
VikS + T]ik 



as follows : 

a l = KjVik — XikVij , 

b\ = XijTjik + Hij^ik — KkVij + ^ikVij , 

(17) c\ = fiijTjik - HikVij 



(18) 



3 For sake of shortness we give the expressions only for some 
coefficients, the other ones being defined in a similar way 
changing j with k or I. 
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and similarly for (02,62^2) (changing k with I) and 
(a 3 ,fe 3 ,c 3 ) (changing X lkl fi ik , v ik , r\ ik with pi, p , g X) 
q ). The system formed by Eqs.(f2l]), @ and (||| is 
closed and may be solved numerically to find the three 
unknown parameters (9 s ,S,9 p ). Then, one of the Eqs. 
@, @, (|l|) and (|0|) may be used to finally get the 
other source coordinate r s . Actually, it is still possible 
to further reduce the system using the resultant method 
( [Erdl fc Schneider, 1993j ) io get only two equations in 
(9 S , Op) which have to be solved numerically. However we 
have preferred to not apply this method since it is difficult 
to implement in a totally automathic code and it leads to 
no significative decrease of the CPU time with respect to 
the strategy of solving directly Eqs. (^I|), (|22] ) and (f23|). 
Beside, some tests have shown us that the two approaches 
lead to the same results. 

Having found the four parameters (r s ,9 s ,S,9 p ), it is 
now straightforward to find the other two lensing potential 
parameters a and [3. To this aim let us consider again 
Eq.(||) written for images i and j; dividing side by side 
one gets : 



3. Solving numerically for potential parameters 



r, COS tf,; - 



1 ' r l-5cos(26 l -28 p y f3 



1 - S cos (29 j - 29 p ) 



(24) 



The same holds for images k and I. Taking the natural 
logarithm of these expressions and solving for a and (3 one 
finally finds : 



a-kibij — aijbki 
bijCki — bkiCij 



1 



Cij(a - 1) 



'>j 



where we have defined : 
-r s cos (9i -9 S ) 



= In 



1 j — r s cos [vj — w s j _ 

1 -5cos(29 t - 29 p ) 
1 - 5 cos (26 j - 29 p ) 



bij = In 



In (n/rj), 



(25) 
(26) 

(27) 

(28) 
(29) 



and similar relation for a k i, bki, cj~i- Note that all the 
quantities entering Eqs.(^5|) and (^6]) are known since im- 
age positions are observed and (r s , 9 S , 8, 6 P ) have been 
found solving the system formed by Eqs.(p7|), (|2l|), ( p2[ ) 
and d23|). These latter equations, Eqs.(p5|) and ( pq ) allow 
us to find out the source positions and the lensing poten- 
tial parameters by knowing only the images positions and 
without using any information on the flux ratios (which 
may be affected by microlensing, that is images cannot be 
resolved) and/or the time delays among images (which are 
measured only for two quadruple systems). 



Eqs. fl2ip, (|22j), (23) are a system of three equations which 
may be solved to find out the source position angle 9 S , 
the flattening indicator 5 and the position angle 9 p of the 
lensing potential. Then Eq. (|17|) gives the other source co- 
ordinate r s and, finally, Eqs.(|25|) and ( p6| ) allow to get the 
other two potential parameters, the slope a of the radial 
profile and the boxiness parameter (3. The system formed 
by the first three equations is not linear and its solutions 
may be found only numerically. This has led us to develop 
an algorithm^ to search for the solutions of the system. To 
avoid introducing any bias in the search, we give to the 
algorithm TV random starting points for (9 S , 5, 9 P ), where 
M is a number fixed by the user^J. We then obtain M so- 
lutions which are not all different from each other and are 
not all physically acceptable. To select among these we 
have imposed a set of selection criteria: the code checks 
the list of solutions and finally retains only the ones 
satisfying the whole set of criteria. Schematically the se- 
lection procedure works as follows : 

1. check that the solutions (9 S , S, 9 p ) are good approxima- 
tions of the real ones, i.e. retain only those solutions 
which inserted again into Eqs.(^l|), ( |22| ) and (|3|) solve 
them with high accuracy; 

2. evaluate r s from Eqs.Q, (jlf), @, © and exclude 
all solutions which give values of r s different among 
each other and/or negative: this test is necessary to 
be sure that the quadruplet (r Sl 9 s ,S,9 p ) indeed solve 
the system (i.e. it is not a false solution due to a bad 
convergence of the numerical code) and to erase the 
unphysical solutions with negative r s ; 

3. select solutions with 5 < 1 : this is necessary to be 
consistent with the starting hypotheses; 

4. look for solutions with a < 2 in order to have S(r, 9) 
monotonically decreasing with the radius; 

5. exclude all solutions with \(3\ > 1/2 since 
t hese often lead to u nphysical mass - radius relation 
( |Zhao fc Pronk, 2001] ); 

6. finally retain only solutions which give rise to galactic 
models with 0.2 < q K < 1.0 where q K is evaluated 
through Eq.(^); this latter constraint has been added 
since galactic halos (which represent more than 90% 
of the lens galaxies mass) are never prolate and their 
flattening (even if not well constrained) is always larger 
than 0.4 so that 0.2 is a quite conservative lower limit. 

Applying this set of constraints eliminates a lot of so- 
lutions leaving us with only A4 solutions, being M. a num- 
ber much lower than M . There are now two possibilities. 

4 The code is nothing more but a notebook written for 
MATHEMATICA, which we have named PULP (Probing 
Unknown Lens Parameters). 

5 N should be large to explore a wide region in the parame- 
ter space, but not too large to save computer time. The right 
choice must be a compromise between these two different cir- 
cumstances. A possible strategy could be to fix N = 4000 and 
then, eventually, run PULP more than one time if necessary. 
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The first is A4 =0; this may be due to a wrong choice 
of the potential, that is the lensing potential does not be- 
long to the class described by Eq.(j(|), but to be sure that 
this is indeed the case, it is better to run the code more 
than one time with M very large in order to explore a 
very wide region of the parameter space. The other hy- 
pothesis is that M. > 1, i.e. we still have more than one 
solution^. After having checked that these M. solutions 
are really different (since it is still possible that one has 
the same solution except for the values of 9 S and/or 9 P 
which, being angles, are defined mod 2tt), we have to fur- 
ther select among these ones. A powerful discriminator is 
the time delay ratio among different couples of images, but 
unfortunately this is not useful in practice since only two 
of the fifteen observed quadruple systems have measured 
time delays^. On the other hand flux ratios are always 
measured with a good accuracy; these quantities depend 
on the lensing potential parameters, so that one may use 
their observed values as constraints to select among the 
M. solutions previously found. 

To this aim, let us first evaluate the flux ratios for the 
class of lensing potentials described by Eq. (||) . Let 4> s be 
the (unknown) source flux and cf>i the one coming from the 
image with coordinates (7^,6^). The magnification due to 
the lensing effect is given by ( |Schncidcr et al., 1992 ) : 



being : 



1 

K = — 

2 



9 2 
dx 2 



7i 



72 



dx 2 



1(1 



9 2 
dy 2 

dy 2 



7? 



7ll 



1 d 2 Tp 

2 dxdy 



(30) 



(31) 



(32) 



(33) 



which are, respectively, the dimcnsionless surface mass 
density, yet introdu ced in Eq.([^), and the t wo components 



of the shear vector ( Schneider et al., 1992 ). 

Introducing Eq.|f) into Eqs.(|3l|), (|32j), (|33j) and then 
the results into Eq.po|), after a lengthy calculation (be- 
ing careful that our transformation to polar coordinates 
is a little bit different from the usual) one finally gets the 
magnification which is 



fi- l (r,e) = |1 - (a 2 + f 2 + f2)rr 2 H&) + 



fi 



(a - 1) (a 2 
being : 

4<5/3[cos (29 - 29 p ) - 6] 



af 2 )r 2{a - 2) F(9)\ 



sm = 



[l- 8 cos {29 -29 



(34) 



(35) 



p) 



6 Unfortunately the case M — 1 is very unlikely. 

7 Detailed informations on the observed lens systems, both 
double and quadr uple, may be found in t he CASTLES collab- 
oration web page (Kochanek et al., 2000). 



The flux ratios between images i and j is then simply : 

(36) 



y 4>i 4>s 4>i v{n,0i) 



which is more convenient to express in magnitudes as : 



= TTtj — mi = —2.5 log< 



'13 ) 



(37) 



being rrii (rrij) the observed magnitudes of image i (j). 

For each solution among the M. previously selected, 
we may now evaluate the flux ratios with respect to the 
image i and compare these with the observed ones finally 
retaining only those solutions such that 



mij(obsd) — m.ij(sol) 



rriij (obsd) 



< e, 



(38) 



being rriij(obsd) the observed value, rriij(sol) the ones 
predicted by the examined solution and e a tolerance fixed 
by the user. The value of e is difficult to choose for two 
main reasons. On the one hand, is extremely sensitive 
to the values of a and f3 since these quantities enter as 
exponents in the theoretical formula. That is why these 
parameters must be finely tuned to be in agreement with 
observations since a small error on their values may lead to 
a high error on rriij. On the other hand, flux ratios may be 
affected by microlensing which may change the flux of only 
one of the images leading to a flux ratio which is different 
from what is predicted using Eqs.(|34|) and (|3^). These 
reasons lead us to not choose a too small value of e in order 
to avoid to exclude solutions not satisfying the constraint 
(|38|), which are indeed good approximations of the right 
one. Some tests have shown that a good choice could be 
e = 0.5 even if a so high value is not too restrictive. Given 
this problems of fine tuning, we have decided to not use the 
flux ratios as constraints when modelling real lenses since 
it is very difficult to reach a so high accuracy in recovering 
the parameters a and (3 when working with data affected 
by observational errors. 

The algorithm searches for solutions and matches the 
whole set of constraints finally leading to a very small 
number of solutions or directly selecting only one solution. 
In order to be sure that this is the correct solution one 
may also recalculate the image positions to see if they are 
in agreement (within the observational errors) with the 
observed ones and finally choose only one solution. 

3.1. A test example 

To test if our code works well or not, i.e. if it is able to 
recover the correct solution solving the system given by 
Eqs.(p"7|), (^l|), ( |22] ) and (p3|), we have made some tests on 
simulated situations. As an example, we describe one of 
these. 

As a first step, we have given the source position 
(r s , 9 S ) as (0.09", 35°) and chosen the lensing potential pa- 
rameters as (a, /3, 5, 9 p ) = (1.3, -0.25, -0.25, 60°.5). Then, 
we have numerically solved lens equations (||) and (|]) for 
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this potential determining image positions. We have also 
evaluated flux ratios with respect to image i and finally 
we have approximated all these quantities to simulate the 
image positions^. These quantities are then used as in- 
put for the algorithm; we have then fixed Af = 1000 and 
then started the search. The first set of selection criteria 
(not considering the one on the flattening) gives 72 (~ 7% 
of the starting ones) solutions, which further reduces to 
24 (~ 2% of the starting ones) after the constrain on q K . 
Finally the conservative constraint on the flux ratios leaves 
us with only one solution : 

(r s ,6 s ,5,9 p ,{3,a) = 

(0.09", 34°.6, -0.24, 240°.5, -0.26, 1.29) 

which is a very good approximation of the starting 
parameters given in the simulation. Note, however, that 
the recovered value of 6 P differs from the real one of ir, 
which is however a well known degeneracy which we do 
not comment on further. 

The test has been repeated many times with different 
sets of source positions and potential parameters and has 
shown that our code always recovers the exact values of 
the parameters with a quite good accuracy. 

4. Application to PG1115+080 and estimate of 
the Hubble constant 

Having checked that the semi - analytical method indeed 
works in recovering the lensing potential parameters, the 
next step is to apply it to real lenses to see whether the the 
results are in agreement with previous ones in literature. 
To this aim, we need a system with four images and a 
good astrometry not only of the image positions, but also 
of the lens galaxy centre since this latter will be used as 
the origin of our coordinate system. Beside, we also need 
that there is only one galaxy acting as lens otherwise it 
is not possible to factorize the lensing potential. In the 
CASTLES database ( |Kochanek et al., 2000] ) there are fif- 
teen quadruple lenses but only some of them satisfy (more 
or less) all our requirements. Among these we have cho- 
sen to consider the well known PG1115+080, first discov- 
ered by Weymann et al. (1980) and then studied in detail 
by several authors (se e, e.g, |Kcctoii fc Kochanck, 1997 : 
Williams fc Saha, 2000| ; fch&o fc Pronk, 2001|) with dif- 
ferent techniques. This system consists of four images 
(named ^41, A2, B and C) of a radio quiet QSO at 
z s = 1.722, while the lens is an elliptical galaxy belong- 
ing to a group of ~ 10 galaxies at zl — 0.310. The cen- 
ter of the group is at (r g ,6 g ) = (20" ± 0.2", -117° ± 3°) 
and its effect has been taken into account as an external 
shear in previous models. This consideration should sug- 
gest that PG1 115+080 is not the best system to apply 



our method since we have neglected any contribution of 
the external shear from the beginning. However, we note 
that previous models never considered boxy potentials^ 
as we do and so it is interesting to explore also this pos- 
sibility. Beside, PG1115+080 is one of the two quadruple 
lenses with measured time delays^] which allows us, on 
the one hand, to use also the time delay ratio to select 
among the possible solutions and, on the other hand, to 
get an estimate of the Hubble constant Hq. In the fol- 
lowing subsections we apply the method to this system 
using as images coordinates the ones measured by Impey 
et al. (1998) with HST observations; the time delay be- 
tween images Al and A2 is too small to be measured, while 
t he one between images B and C is Atsc = 25.0 ± 1.7 d 
( Barkana, 1997 ) with image B arriving last. There are also 
two different measures of the time delay ratio among im- 
ages BC and AB : Schechter et al. (1997) first reported 
tabc = AtAB I Atflc = 0.7 ± 0.3, while a later analysis 
by Barkana (1997) found tabc — 1.13 ±0.18. First, we do 
some general considerations on how to treat the errors. 



4.1. Effect of observational errors 

Our semi - analytical method has been tested on simu- 
lated cases implicitly assuming that the image positions 
are measured with no errors. However, when working with 
real data, one has to take into account also the errors and 
how these affect the determination of the lensing potential 
parameters. Since we do not have anay analytical formula 
which gives the solutions as function of the observed im- 
age coordinates, it is not possible to propagate the errors 
from image positions to lens parameters. However, it is 
still possible to give a qualitative (but quite conservative) 
estimate of these latter quantities. To this aim, we have 
added a 5% uncertainty on the image coordinates of one 
of our simulated systems; then, we have run the code more 
than 100 times each time given randomly the image co- 
ordinates each one within the range obtained adding the 
error. Each one of these run has lead us to an estimate 
of the source coordinates (r s ,9 s ) and of the potential pa- 
rameters (a, P,S,9 P ) which we may compare to the real 
values of the parameters which we have assigned from the 
beginning. These tests have shown us that the code still 
recovers the correct values of the parameters within an er- 
ror which is less than 15%. To be conservative we will add 
a 20% error on the reconstructed parameter when working 
with real systems. 

This error on source coordinates and potential param- 
eters induce an error also on the recalculated image po- 
sitions which we have to estimate. To this aim, we adopt 
a similar procedure solving the lens equations (0) , (^) for 



8 The whole task has been done with a MATHEMATICA 
notebook written by us and named LIGelA (Lensing Images 
Generator Implemented Algorithm) . Both LIGelA and PULP 
are available on request to the authors. 



9 Zhao & Pronk (2001) indeed used also potentials similar 
to the our one, but in their method the boxiness parameter is 
fixed by hand from the beginning. 

10 The other one is B 1608+656 which is too difficult to model 
since there are two lensing galaxies which are probably under- 
going a merging event. 
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Image 


(?", O^j observed 


(?1 ^predicted 


Al 
A2 
B 
C 


(1.173 ± 0.002, 258.78 ±0.10) 
(1.120 ±0.002, 283.05 ±0.10) 
(0.950 ± 0.002, 155.43 ± 0.12) 
(1.397 ±0.002, 40.86 ±0.09) 


(1.19^,250^) 
(1.10lS:S,288l») 
(0.93± ™,l50±i) 
(1.36t° : ?J,43t£) 



Table 1. Observed and predicted quantities for 
PG1115±080. The radial coordinate is in arcsec and the 
angular one is in degrees. 



over 100 realizations of the set of parameters, randomly 
choosing each one in the range (0.8p, 1.2p), being p one of 
the parameters. Then, to be conservative, we give as our 
estimate of the image positions the median of the values 
obtained and as range the ones which contains the 90% 
of the obtained values. Even if not statistically well mo- 
tivated, this procedure should be conservative enough to 
assure us that we are not underestimating the errors. 

4.2. Modelling PG1115+080 and estimates of H 

Here we apply our method to PG1115±080; as first step, 
we estimate the coordinates (r, 9) of the four images from 
the one in Impey et al. (1998) translating the origin of 
the coordinate system from the centre of the C image 
to the galaxy one. A summary of the observed quantities 
is reported in Table 1 where we report also the values 
predicted by our best model. 

This latter has been obtained applying our method 
adopting as input the central values of the range deter- 
mined for each image position. This gives us the solution 
with the following parameters : 

(r„0.) s (0.11", 351°) , 

(a, /3, 5, 8 P ) ~ (1.12, -0.31, -0.31, 325°) . 

To these values we add a 20% error for the reasons ex- 
plained in the previous subsection. Now we have to re- 
calculate the image positions to be sure that this set of 
parameters predicts the correct configuration of the im- 
ages and no other visible images. To this aim we proceed 
as described before generating over 200 set of parameters 
{r s ,6 s ,OL, (3, 8, 9 P ) with each value chosen randomly in a 
range defined as (0.8r s , 1.2r s ) and similar for the other 
parameters. For each of this realizations we compute the 
flattening q K of the surface mass density using Eq.(Q), the 
Hubble constant h and the time delay ratio tabc- These 
two latter quantities may be estimated as follows : 



At 



BC T 100 

2a 



{(2-a)(^-r|) 



- 2(1 - a)r s [r c cos (d c -9 B )-r B cos (9 B - 9 S )}} , (39) 
tabc = {(2 - a)(r 2 c - r|) - 2(1 - a)r s x 



{(2-a)(T^-r|)-2(l-a)r a x 
[r A1 cos (9 A i - 9 S ) - rs cos (9 B - 8 )]}~ 



(40) 



In computing h we adopt a flat cosmology with Q m = 0.3 
and Qa — 0.7 which gives us tioo = 33.37 d arcsec -2 . 
Different cosmological parameters change the estimate of 
h less than 3% so that we can neglect their effect. From the 
list so obtained we erase all the set of parameters which 
generates models with q K outside the range (0.2, 1.0) or a 
value of h outside the range (0.1,1.0) or a negative tabc- 
This leaves us with a consistent set of models each one 
with its own four image positions. The final estimate of 
the coordinates are reported in Table 1 where the central 
value is the median of the distribution and the errors are 
defined such that the range so delineated contains the 90% 
of the values. Similarly we estimate the flattening q K , and 
the Hubble constant which turn out to be : 



a - 52+ 07 
q K — u.oz_ 10 



H a = 56±\l km s -1 Mpc 



-l 



[rc cos ( 



ic 



r B cos iU B 



.)]} 



The mean value of tabc turns out to be 0.79, but the val- 
ues are so sparse that defining a 90% range is meaningless. 
Note, however, that this value is in agreement with the one 
obtained by Schechter et al. (1997). We deserve to conclu- 
sions the discussion of these results. Finally, note that we 
have also found other solutions, but we have deleted them 
since the predicted time delay ratio is grossly inconsistent 
with the estimates of both Schetchter et al. (1997) and 
Barkana (1997). 

5. Discussion and Conclusions 

In this paper we have developed a semi - analytical method 
to reconstruct lensing potential in quadruply imaged sys- 
tems. The method relies only on the image positions and 
allows to recover both the source coordinates and the 
whole set of potential parameters for a broad class of 
non - elliptical boxy potentials. The angular structure of 
the lens potential is a crucial feature of the lens models, 
so it is important to generalize from elliptical to boxy lens 
models to understand how the lens models are affected by 
boxiness. Boxy potentials similar to our models have been 
yet considered by Zhao & Pronk (2001), but our method 
is the first which allows to recover all the model parame- 
ters without giving a priori the boxiness parameter /3. The 
method does not use time delays ratios between any two 
images since these very useful quantities are so difficult 
to measure that nowadays only two quadruple systems 
have been used to measure time delays. The flux ratios 
are used only as broad constraints to avoid the problems 
connected to the possible systematic effect of microlens- 
ing; however, they turn out to be very sensitive to the 
slope a of the radial profile and the boxiness parameter 
(3 so that a fine tuning of these two quantities seems to 
be necessary to fit the flux ratios. The method has been 
tested on simulated cases to see if it works well or not, i.e. 
if it is able to correctly recover the potential parameters. 
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This test has shown that it works indeed being only af- 
fected by the well known degeneracy on the position angle 
Op (i.e. both 6 p and 9 p + it are acceptable solutions). To 
take into account the observational errors in the measure- 
ments of image coordinates, we have used a conservative 
approach to estimate the uncertainties induced on the re- 
covered parameters and how they propagate on the image 
positions predicted by the reconstructed potential. Even 
if not statistically well motivated, our error estimates are 
quite conservative such that we are confident that we have 
not introduced any systematic errors in our reconstruc- 
tion. 

The method described has been applied to the quadru- 
ple lens PG1115+080; for this system we have also used 
the contraints coming from the measured time delay ra- 
tio and estimated the Hubble constant. Our best model 
fits well the four images positions which is a nice result 
especially if one considers that we have not introduced 
any external shear from the group the lens galaxy belongs 
to, as is usually done in the previous analyses of this sys- 
tem. This has allowed us also to recover an estimate of the 
Hubble constant, finding out H = 56+ 12 km s _1 Mpc -1 . 
PG115+080 has been studied in detail by different authors 
using different techniques so that it is interesting to com- 
pare our result with the previous ones in literature. Keeton 
& Kochanek (1997) have used the usual least x 2 paramet- 
ric approach modelling the lensing potential as the sum 
of a term due to the galaxy and an external shear from 
the group. They have considered various elliptical models 
and have found that none of them may fit the image po- 
sitions without any external shear. One of their two best 
fit models describes the lens galaxy as an isothermal (i.e. 
a = 1.0) one with axis ratio q = 0.90 treating the group as 
a point mass. Our best fit model is nearly isothermal (hav- 
ing a = 1.12), but the axis ratio is significantly different 
(being q K = 0.521q;°q). However, a direct comparison of 
the angular profile of the two models is not possible since 
the model by Keeton & Kochanek is an elliptical one so 
that the its boxiness parameter is fixed to be (3 = 0.5, 
whilst our best fit model have (3 — —0.31 which is quite 
different. Actually, Keeton & Kochanek start from lens 
mass models which are suggested by the light profile of 
the galaxy. Our approach aims at reconstructing the lens- 
ing potential; the result cannot be immediately compared 
to the light profile since the potential may be also gen- 
erated by a dark halo embedding the visible component 
of the lens galaxy. It is thus difficult to say if the boxi- 
ness of the potential and the axis ratio of the surface mass 
density refer to the dark or visible component which lead 
us to not state any definitive conclusion from the com- 
parison with Keeton & Kochanek model. It is however 
interesting to compare the results on the Hubble constant 
with the previous ones in literature since a net discrepancy 
could suggest that something is wrong with our method. 
A carefull analysis of the possible systematic effects and 
of the different models suited for PG1 115+080 finally lead 
Keeton & Kochanek to estimate the Hubble constant as 
H a = 51^13 km s -1 Mpc" 1 . This value is in quite good 



agreement with the our one, H = 56+ 12 km s _1 Mpc -1 , 
which is a very encouraging result since we have used no 
contribution from external shear. Zhao & Pronk (2001) 
have fitted a similar boxy model to describe PG1115+080, 
but a direct comparison with their result is not possible 
since they do not quote neither a nor 5, but only what 
they called the effective slope of the radial profile and a 
typical flattening. Beside, they fix a priori the boxiness 
parameter, while in our method this is a reconstructed 
parameter. The disagreement between our and their val- 
ues of the source coordinates may so be a consequence 
of the different boxiness parameters; however, we note 
that we predict a value for H which is in agreement with 
what is reported in the second column of their Table 2. 
Finally, it is interesting to compare our result to the one 
obtained by Williams & Saha (2000) using the fully non 
parametric pixelated lens method. In this approach, one 
does not try to reconstruct the lensing potential, but only 
the time arrival surface which, as they show, it is enough 
to estimate the Hubble constant. Combining the results 
from both PG1115+080 and B1608+656 they finally quote 
H = 61 ± 18 km s -1 Mpc -1 (90% confidence limit), still 
in agreement with what we have obtained. 

All these results are encouraging and seem to sug- 
gest that boxy potentials may describe quadruple lenses 
without the need of introducing an external shear. It is 
well known that this latter term is often necessary to fit 
well the image positions in quadruple systems and, be- 
side, Witt (1996) has analytically shown that the min- 
imum amount of shear, necessary to fit quadruple sys- 
tems, is not trivial if the lensing potential is elliptical. 
However, we stress again that our models are not ellip- 
tical so that the result obtained by Witt may not be ap- 
plied. Actually, we succeed in fitting the image positions in 
PG1 115+080 also without introducing the external shear 
and have also predicted a value for the Hubble constant 
which is in good agreement with the previous estimates in 
literature obtained from this lens with different methods. 
These results are encouraging and we are now carrying on 
a detailed study of other quadruple lenses to see whether 
boxy potentials may fit the data and whether the pre- 
dicted galaxy models are physically motivated (Cardone 
et al., 2001). On the other hand, it could be also necessary 
to consider boxy potentials with the external shear since 
many lens galaxies belong to a cluster whose effect should 
be considered before drawing conclusive results. This con- 
sideration lead us to further improve our method to take 
into account also the external shear. This task is not sim- 
ple since to add the shear to the potential introduces other 
two variables, the shear strength 7 and its position angle 
Oj, so that we will have eight variables. Image positions 
give us eight equations so that one may still hope to solve 
the problem since it is still possible to write a closed sys- 
tem of equations which we are tryng to manage to reduce 
to a numerically solvable one. 
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